Under consideration for publication in J. Fluid Mech. 119 

Entropic Lattice Boltzmann Method for 
Large Scale Turbulence Simulation 

By I. V. KARLIN 1 f S. ANSUMALI 1 E. DE ANGELIS 1 J 

H. C. OTTINGER 1 

AND 

S. SUCCI 2 

1 ETH-Ziirich, Department of Materials, Institute of Polymers, ETH-Zentrum, Sonneggstr. 3, 

ML J 19, CH-8092 Zurich, Switzerland 

2 Istituto Applicazioni Calcolo, CNR, viale Policnico 137, 00161 Roma, Italy 

(Received 30th May 2003) 

Recently, a minimal kinetic model for fluid flow, known as entropic lattice Boltzmann 
method, has been proposed for the simulation of isothermal hydrodynamic flows. At vari- 
ance with previous Lattice Boltzmann methods, the entropic version permits to describe 
the non-linear dynamics of short scales in a controlled and stable way. In this paper, we 
provide the first numerical evidence that the entropic lattice Boltzmann scheme provides 
a quantitatively correct description of the large-scale structures of two-dimensional tur- 
bulence, free of numerical instabilities. This indicates that the entropic lattice Boltzmann 
method might provide a starting basis for the formulation of a new class of turbulence 
models based on genuinely kinetic principles. 



1. Introduction 

{} The dynamics of fully developed turbulent flows is characterized by the simultaneous 
non- linear interaction of a large number of degrees of freedom, far too many for analytical 
treatment and also for the processing capabilities of even most powerful computers. 
Since most turbulent flows of practical relevance cannot be simulated ab initio, they 
must be modeled. The goal of such a modeling is to predict the effects of unresolved 
scales of motion, those that cannot be represented within the computer simulation, on 
the resolved ones. Mathematically, these effects are represented by the divergence of the 
Reynolds stress tensor, T — (u'u'), where u' denotes the short-scale component of fluid 
motion, and the pointed bracket stands for ensemble averaging. 

A popular class of models is based on the notion of 'eddy-viscosity', according to which 
the effects of short scales can be incorporated within an effective viscosity v c = v + Vt 
acting upon the large scales. Here v is the bare, 'molecular' viscosity, whereas vt denotes 
the turbulent viscosity. The effective viscosity v e is then expressed as a function(al) of 
resolved fields. The simplest option consists of choosing a suitable algebraic relation 
v t — v t (S), where S is the magnitude of the shear tensor S = VU + V(U) T , U denoting 
the resolved flow field. The simplest eddy- viscosity models express the Reynolds stress 
tensor in an algebraic form: T = i/ e (S)S. This formulation has the significant advan- 
tage of leaving the form of the Navier-Stokes equation for the large eddies intact, with 
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an effective viscosity v e replacing the molecular viscosity v. However, such an algebraic 
representation often fails to reproduce strongly off-equilibrium turbulent effects, such as 
those occurring in the vicinity of solid walls. A better account of these phenomena is ob- 
tained by linking v e to the actual content of turbulent kinetic energy k =< u' 2 /2 > and 
energy dissipation rate e = dk/dt, via v e ~ k 2 /e . Even so, important phenomena such 
as turbulent gusts triggered by local i nstabilities (corresponding to negative effective vis- 
cosity) cannot be captured l)Kraichnanlll976|) because u e is bound to be positive definite. 
Modern formulations which allow for negative viscosities and corres ponding backscatter 
effect s in fc-space have indeed been developed in the last two decades ijLesieur and Metaisl 
However, to the best of our knowledge, none of these models has gained universal 
consensus to date. 

All eddy-viscosity models rely on the analogy between the turbulent transport and 
the molecular transport. Within this analogy, the small-scale eddies play the role of the 
molecules, while the smallest coherence length, Ik, known as Kolmogorov length, plays 
the role of the mean-free path. It is then natural to define a 'turbulent Knudsen number' 
for an eddy of size I as 

Kn t (l) = l j 

It is then understood that small-scale eddies of size close to If. are in local equilibrium 
with eddies of size I, so long as these two scales are well separated, i.e. Kn t {l) « 1. 

A major problem with this analogy is that turbulence features a continuum spectrum 
of excitations, so that the above scale-separation argument simply does not hold. In 
particular, from the above definition of turbulent Knudsen number, it is clear that the 
local equilibrium assumption becomes less and less tenable as I — > Ik- This lack of scale 
separation is a fundamental problem which lies the heart of all failures to develop a 
consistent theory of fully developed turbulence. 

From the computational point of view, the Kolmogorov length is replaced by the grid 
spacing 5 >> Ik, indicating that the dynamics of the smallest resolved scales faces a 
situation similar to finite-Knudscn flows at Kn ~ 1. 

It has been recently pointed out that the kinetic representation of hydrodynamics 
provides a natural generalization of the no tion of eddy- viscosity to such non-equilibrium 
high- if n t regimes ijChen et alJll999l 12003). The key point is that solutions of the kinetic 
equation apply at all orders of the Knudsen number, so that any kinetic model ensuring 
correct hydrodynamic behaviour would handle the dynamics of small eddies at Kn t ~ 1 
in a way which goes beyond the eddy-viscosity representation. 

The practical problem with kinetic theory, however, is the enormous increase of de- 
grees of freedom associated with the (true) Boltzmann equation. Thus, in order to apply 
the aforementioned kinetic approach to fluid turbulence, drastically simplified versions 
of the Boltz mann equat i on need to be dev eloped. The lattice Boltzmann method (here- 
after LBM) Ipucc i 2001; B enzi is a good candidate for such a task. Indeed, 
the method has been applied to a large var iety of fluid flows, including turbulent ones 
l|Benzi and Sucdlll99d: iMartinez efalM mA) . However, the standard lattice Boltzmann 
method meets with severe difficulties in handling the dynamics of near-grid scales with 
I ~ S. In particular, owing to the lack of a H-theorem, these scales often exhibit disruptive 
non- linear instabilities. 

Recently, a new approac h to stabilize the kinetic scheme based on thermo d ynamic con- 

sider ations was developed ijKarlin et a/.ll999llAnsumali and Karlinl2002albl lci lBoghosian et al\ 
l200l|) . The modified method known as "entropic Lattice Boltzmann method" might add 
a further boost towards the formulation of a kinetic model of fluid turbulence. In par- 
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ticular, its non-linear stability, and the ensuing positivity of the distribution function, 
automatically enforces an important realizability constraint which is missing in stan- 
dard Lattice Boltzmann representations. More specifically, this model provides a local, 
adaptive, regulator of the relaxation time, which can be regarded to all effects and pur- 
poses as a turbulence model inspired by genuinely kinetic requirements. In this paper, 
we present a set of numerical experiments which are intended to test this conjecture in 
a semi-quantitative sense. 

The work is organized as follows: first, a description of the LBM and its entropic 
version is introduced Sec.|2J). Then, simulations of two dimensional decaying turbulence 
are presented for the fully resolved as well as unresolved case (Sec. (3J|. For validation 
purposes, results are compared with spectral simulation of the Navier-Stokes equations 
in the fully resolved case. 



2. Lattice kinetic theory 

In the last decade, it has become apparent that minimal versions of the Boltzmann 
kinetic equation provide a fairly efficient alternative to the discretisation of the Navier- 
Stokes equat ions for the simulation of a variety of complex flows, including turbulent ones 
ljSiiccill200lf) . In its simplest instance, this minimal kinetic equation takes the following 
form: 

d t fi + Ci ■ d^h = -T- 1 (Ji - .O , (2.1) 

where /,; = f(x,Ci,t) denotes the probability of finding a particle at position x and time 
t, moving along the discrete direction Cj. The right-hand-side of this equation represents 
collisional relaxation to a local equilibrium on a time scale r. The discrete velocities 
must be chosen in such a way as to ensure sufficient symmetry to recover the basic mass, 
momentum and momentum-flux conservations, so that the Navier-Stokes equations are 
recovered as the large-scale limit of the discrete kinetic equation. For this purpose, many 
options have become available over the years. The major appeal of the kinetic approach is 
its remarkable simplicity, which translates into a corresponding computational efficiency. 
A distinctive mark of this simplicity is the fact that, i) the streaming operator is linear 
and proceeds along a set of constant directions, ii) the non-linearity is entirely due to 
the local equilibrium, which is fully local in configuration space. This is in vivid contrast 
to the Navier-Stokes representation in which both non-locality and non-linearity are 
condensed into a single term, namely the convective operator u ■ 9 x u. Typically, the local 
equilibrium is chosen in the form of a quadratic polynomial in the fluid speed : 



ft q = P w t 



kiCi ■ u + fc 2 (c, ■ u) 2 + fc 3 (ciCi) : (uu) 



(2.2) 



where, Wj, fci, k2, and are lattice dependent constants depending on the lattice sound 
speed c s (a constant for the present case of athermal flows). The advantage of this ap- 
proach is that local equilibria can be readily constructed after particles have gone through 
the streaming phase. Subsequent irreversible relaxation to this local equilibrium provides 
viscous behavior, with a kinematic viscosity of the order of v ~ c 2 s t. The simplicity of 
this scheme is hard to beat. However, the fact that local equilibria are prescribed at 
the outset and do not result from a self-consistent relaxation dynamics in kinetic space, 
implies that compliance with the second law of th ermodynamics, the H-theorem in ki- 
netic language, is generally lost l|Succi et aZJ l2002f) . The result is that, whenever large 
gradients develop on the lattice scale, as it is the case for fully developed turbulence, 
standard LBE schemes are subject to numerical instabilities. This provides a strong mo- 
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tivation to develop a new class of lattice Boltzmann schemes, capable of accommodating 
the H theorem. A particularly interesting option has re cently emerged in the form of the 
so-called entropic Lattice Boltzmann method (ELBM) ijAnsumali and Karlinll2002cj) . In 
remainder of this section we will describe entropic model. 

The basic strategy of the entropic Lat tice Boltzmann method is to write the dynamic s 
in terms of a properly chosen H function l)Karlin et oill999llAnsumali and Karlinl2 002c'l . 
For isothermal hydrodynamics, the discrete H function is found to be in Boltzmann-form: 



H {w i ,c i} = Y,fi^[£r), (2.3) 



A 

where b D is the number of discrete velocities in ^-dimensions. The discrete velocities are 
tensor products of the discrete velocities in one dimension and the weights are constructed 
by multiplying weights associated with each direction. The minimal set of discrete veloci- 
ties needed to reconstruct the Navier-Stokes equations are related to zeroes of third-order 
Hermite polynomials in a and for D = 1, the three dis crete velocities are c ,; = c{ — 1, 0, 1}, 
whereas the corresponding weights are w = {i, §, ^j llShan and Helll99^ . The discrete- 
velocity local equilibrium is the minimizer of the corresponding entropy function under 

the constraint of local conservation laws: 53j=i /f q {L c i} = {p, P u }- The explicit solu- 
tion for the /j° q is 



D 



pWiJl 



2 _ l/1+ ^)to^f\ ( , 4) 



with j being the index for spatial directions. Further, a notion of the bare collision A, 
defined as the collision term stripped of its relaxation parameters is introduced. In the 
case of BGK model A = f oq — f , namely the bare departure from local equilibrium. The 
time stepping in this method is performed through an over-relaxation collisional process 
and linear convection through a sequence of steps in which the H function is bound 
not to decrease. The monotonicity constraint on the H function is imposed through the 
following geometric procedure: In the first step, populations are changed in the direction 
of the bare collision in such a way that the H function remains constant. In the second 
step, dissipation is introduced and the magnitude of the H function decreases. Thus, 



f i (-x,5t) = f i (x-c i St,O) + a0 



(2.5) 



were (3 is the discrete form of the relaxation frequency related to r 1 and the parameter 
a is defined by the condition: 

H(f) = H(f + oA), (2.6) 

and close to the local equilibrium a is equal to 2 ijAnsumali and Karlinll2002a(l . The 
scheme is illustrated graphically in figure ^ As shown in figure ^ in order to find the 
desired point f (/?), we first find the point f* on the curve L of constant H function. This 
way of enforcing the H theorem ensures non-linear stability. A n algorithm to implement 
Eq. (|2.6|) has been presented bv lAnsumali and KarlirJ j200 2a). The local adjustments of 
the relaxation time (via the parameter a), as dictated by compliance with the H-theorem, 
guarantee positivity of the distribution function also for the case of discrete steps, thereby 
ensuring non-linear stability of the numerical scheme. 

A comment on the relation between non-linear stability and built-in subgrid modelling 
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Figure 1. Stabilization procedure. Curves represent entropy levels, surrounding the local equi- 
librium f oq . The solid curve L is the entropy level with the value H(f) = H(f*), where f is 
the initial, and f* is the auxiliary population. The vector A represents the collision integral, 
the sharp angle between A and the vector — V-ff reflects the entropy production inequality. 
The point f* is the solution to Eq. 12.61 . The result of the collision update is represented by 
point f(/3). The choice of /3 shown corresponds to the 'over-relaxation': H(f(/3)) < H(f) but 
H(f(/3)) > H(M). The particular case of the BGK collision (not shown) would be represented 
by a vector Abgk, pointing from f towards f eq , in which case M = f cq . 



is in order. As is well known, in a turbulent flow energy (enstrophy in two-dimensional 
set up) cascades from large to small scales. While physically the cascade picture is ter- 
minated by dissipation at the Kolmogorov length Ik , in an underresolved simulation, the 
cascade needs to be truncated at the grid spacing 5 > Ik- As the entropy of the ELBM 
("grid entropy") accounts for the net effect of all degrees of freedom that are not re- 
solved in the simulation (including the physical entropy), it provides the natural key to 
implementing artificial dissipation on the grid scale. According to the H theorem, the 
cascade is terminated on the grid scale in such a way that the subgrid scales cannot "fire 
back" at scales larger than 8. Consequently, the grid entropy of the ELBM allows us 
to decouple all subgrid effects from the cascade in a most natural way. This is exactly 
the goal traditionally approached with the concept of eddy-viscosity by pushing subgrid 
effects to the decoupled level of local equilibrium thermodynamics. 



3. Decaying turbulence 

In this section, we compare the results of the Entropic Lattice Boltzmann simulation 
with those of pseudo-spectral simulation of the Navier-Stokes equations (hereafter SS) 
for the case of two-dimensional, homogeneus, incompressible decaying turbulence. In this 
setup, one starts with a Gaussian random field in a periodic box, and probes the decay 
of the turbulent field. 

The initial conditions are given by a zero-mean Gaussian vorticity field with random 
Fourier phases. The functional form of the energy spectrum is frequently chosen as: 

E(k)=C k A [l + (k/k ) B+1 r\ (3.1) 

where, Co is a normalization constant and parameters A, B and fcn are cho sen in such a 
way that the energy spectrum remains narrow banded i Bracco et qjJl2f)f)(T^ . 



As a SS algorithm, we use the pseudo-spectral method, i.e., the non- linear terms are 
evaluated in real space, with a formulation in primitive variables (the two velocity com- 
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ponents). The equation to be solved in Fourier space looks as follows: 

du ' " [ 

-— ^ + jkj(ujUi) = -jhp + -^k 2 Ui (3.2) 

where the symbol Q denotes the Fourier transform of the quantity Q. At each time 
step the nonlinear term is projected onto a solenoidal field in order to maintain in- 
compressibility condition kiiii = 0; this procedure is straightforward in /c-space, being 
V(jkj(ujUi)) = jkj (8u + UiUi/k 2 ) (ujUi), where V is the projecting operator. The eval- 
uation of the nonlinear terms is performed with a 3/2-dealiased pseudo-spectral method. 
For the time-mar ching, a, semi -implicit low storage third-order accurate Runga-Kutta 
scheme was used arn 11988). 

In the following subsections, a detailed comparison of the result of decaying turbulence 
simulated with SS and ELBM is presented. 

3.1. Fully-resolved simulation 

We consider a two-dimensional box of length L = 512 and viscosity v = 5.0 x 10~ 4 (all 
quantities are in lattice units), for the SS as well as the ELBM. As described earlier, 
the initial conditions are given in Fourier space with ko = 0.2, A = 6 and B = 17 (see 
Eg. 13. 1|) . The initial spectrum is reasonably peaked at low wave numbers. The initial 
velocity profile has mean kinetic energy E(t = 0) = 1.64505947 x 10 -4 and the mean 
enstrophy is Z(t = 0) = 2.0348495 x 1 0~ 6 . A rough estimate of the eddy turnover time is 
t e w Z- 1 ' 2 w 700 jBracco et a/.ll2000|) . The Reynolds number based on the mean initial 
kinetic energy and the box- length as the characteristic length is Re = L-J (2E)/u w 
13134. The estimated dissipation length is Ld ~ LRe~ x l 2 w 3.2. 

The agreement between the two simulations is very good up to t/t e ~ 20. As, shown in 
figure |21 at t/t e ~ 30 ( i.e., t ~ 20000) slight difference between the two methods starts 
to show up in the vorticity contour plot. However, this difference is only in the small 
scale features of the flow, while all large scales are still the same at this time. 

Figure [21 shows a comparison of the energy and enstrophy spectra, respectively, up 
to t = 10000. The plot shows that the energy and the enstrophy are narrow banded 
initially. The enstrophy and the energy plots for the SS show the typical behavior, where 
during the course of simulation high wave number modes start gaining enstrophy (and 
consequently some energy too). 

In freely decaying 2D turbulence, dissipation takes place via a enstrophy cascade from 
large to small scales. To prevent enstrophy pile-up, a small-scale enstrophy sink is typi- 
cally required in numerical simulations. Figure|3shows that the ELBM naturally provides 
a cutoff at high wave-number. In the subsequent section, we shall show that this cutoff 
does not affect the low wave- numbers region in any unphysical manner. 

Figure 01 shows the energy and enstrophy time decay respectively. We see that these 
quantities are slightly underestimated by the ELBM. This is probably related to the 
sharp cutoff at high wavenumber present in the ELBM simulation. 

As the simulation is performed in the low-Reynolds number regime, both methods 
provide a rather poor agreement with the Kolmogorov-Batchelor theory, and show visible 
deviations from the theoretically expected spectrum E(k) ~ fc~ 3 . 

3.2. Unresolved simulation: example I 

In the previous section, 13.11 we have verified that the entropic method can be used to 
perform a fully resolved simulation of turbulence. We also observed that when the spec- 
trum is not fully resolved by the computational domain, the ELBM method introduces a 
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Figure 2. Contour plots of vorticity at t = 0, 5000, 20000: ELBM results (right column), 
spectral results (left column). Approximate eddy turnover time is t e = 700. 




k k 

Figure 3. Energy and enstrophy spectra (left and right plot respectively) in the fully resolved 
simulation (example II) for three different times: t = 1000, 5000, 10000. Solid line is for spectral 
simulation and the dotted symbols represent ELB simulation. A dashed line showing k spectra 
in the energy spectra is also shown. 

sharp dissipation at high wave numbers. In this section, we shall explore how such high-/c 
cut-off affects the spectrum of the resolved scales. 

In order to magnify this effect, for the spectral simulation we used 256 x 256 grid points, 
while for the ELB we performed the simulation with 128 x 128 grid points. We prepare 
the initial condition in such a way that the resolved scales are the same for both methods 
(see plot on the left hand side in figure EJ- The spectral simulation is fully resolved while 
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Figure 4. Time history of the energy (left plot), and the enstrophy. Eq and Zo are the energy 
and the enstrophy respectively at t — 0. Solid line SS and dotted line ELBM. 




k k 

Figure 5. Wave number energy spectra in the unresolved simulation (example I). The plot 
on the left hand side shows the initial condition, while that on the right hand side shows the 
energy spectra for three different times: t = 4000, 20000, 40000. The solid line is for the spectral 
simulation and the dashed symbols represent ELB simulation 

the spectrum for the ELBM simulation is truncated almost in the middle. The plot of 
the energy spectra., figure 5, at. different times shows that the ELBM simulation follows 
the spectral simulation reasonably well. At low wave- numbers, i.e., in the resolved scales, 
the deviation from the spectral simulation is quite small. Near the cutoff wavenumber, 
a slight (but localized) hump in the profile obtained from ELBM appears. This hump 
and the agreement at low wavenumbers, can be interpreted as a signature of the very 
localized nature of the ELBM dissipative cut-off at high wavenumbers. 

3.3. Unresolved simulation: example II 

In the previous two sections, we showed that the entropic method can perform turbulence 
simulation, both resolved and underresolved, to a reasonable degree of accuracy. In this 
section, we shall explore the ELBM behaviour for flows much beyond the grid resolution. 
For this case, we have chosen the box-size as L = 512 and viscosity v = 5.0 x 10~ 6 for 
the ELBM simulation. As described earlier, the initial conditions are given in Fourier 
space with fco = 0.2, A = 6 and B = 17 (see Eg. 13. If) . The initial velocity profile has 
mean kinetic energy E(t = 0) = 1.20861512 x 10~ 4 and mean enstrophy as Z(t = 0) = 
1.4949915 x 10~ 6 . A rough estimate of the eddy turnover time is t e w Z~ 1/2 w 818 
l|Bracco et alj f2000). The Reynolds number based on the mean initial kinetic energy 
and the box-length as the characteristic length is Re = Ly/(2E)/v « 1.59 x 10 6 . The 
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Figure 6. Energy and enstrophy spectra (left and right plot respectively) in the ELBM unre- 
solved simulation (example II) for three different times: t = 5000, 15000, 20000. A dashed line 
showing k~ 3 spectra in the energy spectra is also shown. 

dissipation length scale is estimated as L^ ~ L Re^ 1 ^ 2 ^ w 0.4. Figure[f)]shows that during 
the time evolution, Batchelor spectra E). ~ fc~ 3 for the energy and for the enstrophy 
in the inertial regime, are reproduced to a reasonable accuracy. The scaling law is best 
fitted at t = 15000 as 2 x 10~ 7 fc~ 3 . A rough estimate of the enstrophy transfer rate 
yields 1.7 x 10~ n . This gives a value of 4.8 for the constant appearing in the Batchelor 
scaling, while the theoretical value is 2.626. We remind that the determination of this 
constant requires an accurate evaluation of the enstrophy transfer rate, so that we expect 
the agreement to improve with increasing system size. The agreement with the scaling 
law in such an underresolved simulation shows that the ELBM exhibits a sort of built-in 
subgrid model for turbulent flows. 

It is argued that in two-dimensional isotropic turbulence, the eddy- viscosity develops a 
dependence on the wavenumber and m ay depart from its mean value in both positive and 
negative directions ijKraichnanl 1 1 976| ) . It is therefore of interest to investigate whether 
ELBM has similar features. We can argue that the deviation of the effective relaxation 
frequency a/3 from its near equilibrium value 2/3 is related to the deviation of viscosity 
from the molecular viscosity v{0) — (1/2/3 — 1/2). In order to parameterise the relax- 
ation behavior, we define the 'effective viscosity' of the ELBM fluid (in lattice units) as: 



In figure {7\ the probability distribution of turbulent fluctuations in effective viscosity as 
compared to the bare molecular viscosity, bvjv = v c g(a, (3) / V(/3) — 1 is plotted. From 
this figure a three-modal distribution is clearly visible, corresponding to three distinct 
classes of events: 

1) Plus-events, [Svjv > 0): Effective viscosity is much larger than molecular viscosity. 
2) Normal-events, (5v/v = 0): Effective viscosity is identical to molecular viscosity, within 
the machine accuracy. 3) Minus-events, (5v/v < 0): Effective viscosity is much smaller 
than molecular viscosity. This case corresponds to local micro-instability. 

The relative weights, {Mpi us , MNormai, A^Minus}, for the three events at three different 
times, t = 1000, 5000, 10000, are {0.227002, 0.754618, 0.018380}, {0.099759, 0.853831, 0.046590}, 
and {0.097823,0.839673,0.062504}, respectively. This shows that, the number of Minus- 
events grows in time, while Plus-events do the opposite. The above figures have been 
obtained by defining Normal-events via the condition \5v/v\ < 0.01. Such a low thresh- 
old indicates that the distribution of Normal-events is a Dirac's delta. This shows that 




(3.3) 
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Figure 7. Distribution M of 8v/v at time t = 1000, 5000, 10000 in a unresolved simulation. 
Notice the delta function at Su/u = 0, where the distribution is normalised in such a way that 
area under curve is unity. 

the parameter a is highly 'intermittent', e.g. quiescent in most part of the flow, with 
strong bursts of activity in small regions (this is confirmed by visual inspection of the 
spatial distribution a(x, y), not shown here). This type of behaviour cannot be captured 
by any 'conventional' picture of eddy-viscosity. 

4. Conclusions 

In conclusion, we have shown that entropic lattice Boltzmann models (ELBM) can cope 
with the instabilities which plague standard lattice Boltzmann schemes in a regime where 
subgrid scales are dynamically excited. This boost of numerical stability is the result of 
securing compliance with the H-thcorcm, hence positivity of the discrete distribution 
function. Therefore it can be concluded that ELBM exhibit a built-in subgrid model 
which is rooted in a basic principle of statistical mechanics rather than on heuristic 
requirements or high-order discretizations of the Navier-Stokes equations. This built- 
in model cannot be reduced to an algebraic relation between eddy-viscosity and the 
large scale shear. This is not surprising, since the present ELBM accounts for both 
positive and negative turbulent relaxation, the latter corresponding to local instabilities 
which escape traditional eddy-viscosity turbulence models. Future work will tell whether 
such encouraging qualitative features can be turned into a quantitative genuinely kinetic 
approach to the problem of turbulence modeling. 
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